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Abstract 

This work addresses the problem of optimal pricing and hedging of a European 
option on an illiquid asset Z using two proxies: a liquid asset S and a liquid European 
option on another liquid asset Y. We assume that the S-hedge is dynamic while the 
Y-hedge is static. Using the indifference pricing approach we derive a HJB equation 
for the value function, and solve it analytically (in quadratures) using an asymptotic 
expansion around the limit of the perfect correlation between assets Y and Z. While 
in this paper we apply our framework to an incomplete market version of the credit- 
equity Merton's model, the same approach can be used for other asset classes (equity, 
commodity, FX, etc.), e.g. for pricing and hedging options with illiquid strikes or 
illiquid exotic options. 



1 Introduction 

Consider a trader who wants to buy or sell a European option Cz on asset Z with maturity 
T and payoff Gz- The trader wants to hedge this position, but the underlying asset Z is 
illiquid. However, some liquid proxies of Z are available in the marketplace. First, there 
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is a financial index (or simply an index) S (such as e.g. S&P500 or CDX.NA)Q whose 
market price is correlated with Z. In addition, there is another correlated asset Y which 
has a liquidly traded option Cy with a payoff Cy similar to that of Cz, and with the 
same maturity T. The market price py of Cy is also known. 

Our trader realizes that hedging Z-derivative with the index S alone may not be 
sufficient for a number of reasons. First, she might be faced with a situation where 
correlation coefficients p yz ,p sz (which for simplicity are assumed to be constant) are such 
that p yz > p sz . In this case we would intuitively expect a better hedge produced by using 
Y or Cy as the hedging instruments. Second, if we bear in mind a stochastic volatility-type 
dynamics for Z, the stochastic volatility process may be "unspanned", i.e. the volatility 
risk of the option may not be traded away by hedging in option's underlying^} If that 
is the case, one might want to hedge the unspanned stochastic volatility by trading in a 
"similar" option with on the proxy asset Y. So our trader is contemplating a hedging 
strategy that would use both S and Y. To capture an "unspanned" stochastic volatility, 
the trader wants to use a derivative Cy written on Y rather than asset Y directly. 

As transaction costs are usually substantially higher for options than for underlyings, 
our trader sets up a static hedge in Cy and a dynamic hedge in St- The static hedging 
strategy amounts to selling a units of Cy options at time t = 0. An optimal hedging 
strategy would be composed of a pair (a*, it*) where a* is the optimal static hedge, and 
7r* (where < s < T) is an optimal dynamic hedging strategy in index S t . The pair 
(a*, should be obtained using a proper model. The same model should produce the 
highest/lowest price for which the trader should agree to buy or sell the Z-option. 

In this paper we develop a model that formalizes the above scenario by supplementing 
it with the specific dynamics for asset prices St, Y t and Z t , and providing criteria of 
optimality for pricing options Cz- For the former, we use a standard correlated log- 
normal dynamics. For the latter, we employ the utility indifference framework with an 
exponential utility, pioneered by Davis [3J, Hodges & Neuberger [7j and others, see e.g. 
[6] for a review. As will be shown below, this results in a tractable formulations with 
analytical (in quadratures) expressions for optimal hedges and option prices. 

As the above setting of pricing and hedging an illiquid option position using a pair 
of liquid proxies (e.g. a stock and an option on a different underlying) is quite general, 
one could visualize its potential applications for various asset classes such as equities, 
commodities, FX etc. For definiteness, in this paper we concentrate on a problem of 
practical interest for counterparty credit risk management [^J Namely, we consider the 

■"■Here we refer to this instrument as an index, but it could be any "linear" instrument such as stock, 
forward, etc. 

2 For a discussion of such scenarios for commodities markets, see [19] . 

3 Most of the formulae below, excluding those that use specific forms of payoffs, are general and 
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problem of pricing and hedging an exposure to a counterparty with an illiquid debt, and 
in the absence of liquidly traded CDS referencing this counterparty. For such situation, 
no market-implied spreads are available for the counterparty in question. Instead, one 
should rely on a model to come up with theoretical credit spreads for the counterparty. To 
this end, we use a version of the classical Merton equity-credit model [H] which is set up 
in a multi-name setting, and under the physical (i.e. "real", not " risk- neutral" ) measure. 
Most importantly, unlike the classical Merton's model, we do not intend to use firm's 
equity to hedge firm's debt. Instead, illiquid debt is hedged with a proxy liquid debt, 
and a proxy credit index. In what follows, to differentiate our framework from that of 
the classical Merton model, we will refer to it as the Hedged Incomplete-market Merton's 
Dynamics, or HIMD for short. 

1.1 Relation to previous literature 

Our model unifies three strands in the literature on indifference pricing. 

The fist strand deals with hedging an option with a proxy asset, as developed in Davis 
[I], Henderson & Hobson [5], Musiela & Zariphopoulou [15], and others. In this setting, 
one typically hedges an option on an illiquid underlying with a liquid proxy asset. 

The second strand develops generalizations of the classical Merton credit-equity model 
to an incomplete market setting. Typically, this is achieved by de-correlating asset value 
and equity price at the level of a single firm, see e.g. Jaimungal & Sigloch [TT], T. Leung 
& Zariphopoulou [T5]. As long as we do not use firms equity to hedge firm's debt but 
instead use a liquid proxy bond as a hedge, such modification of the Merton model is not 
needed in our setting. 

The third strand is presented by [§j who develop a static-dynamic indifference hedging 
approach for barrier options. Ilhan and Sircar considers hedging a barrier option under 
stochastic volatility using static hedges in vanilla options on the same underlying plus 
a dynamic hedge in the underlying. This results in a two-dimensional Hamilton- Jacob i- 
Bellmann (HJB) equation. Our construction is similar but our hedges are a proxy asset 
and a proxy option, while volatility is taken constant for simplicity. 

2 Static hedging in indifference pricing framework 

Borrowing from an approach of [8] for a similar (but not identical) setting, we now show 
how the method of indifference utility pricing can be generalized to incorporate our sce- 
nario of a mixed dynamic-static hedge. 

applicable for other similar settings. 
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To this end, let H(Y T , Z T ) be the final payoff of the portfolio consisting of our option 
positions, i.e. 

U a (Y T ,Z T ) = G z -a*G Y (1) 

As long as both European options Cz,Cy pay at the same maturity T, we can view this 
as the payoff of a combined ("static hedge portfolio") option g%, which involves payoffs 
Gz and Gy of both derivatives Cz and Cy . Such option may be priced using the standard 
utility indifference principle. The latter states that the derivative price is such that 
the investor should be indifferent to the choice between two investment strategies. With 
the first strategy, the investor adds the derivatives to her portfolio of bonds and stocks 
(or indices^]) S, thus taking from, and adding apy to her initial cash x. With the 
second strategy, the investor stays with the optimal portfolio containing bonds and the 
stocks/indices. 

The value of each investment is measured in terms of the value function defined as 
the conditional expectation of utility U(Wt) of the terminal wealth Wt optimized over 
trading strategies. In this work, we use an exponential utility function 

U(W) = -e~^ w (2) 

where 7 is a risk-aversion parameter. In our case, the terminal wealth is given by the 
following expression: 

W T = X T + U a (Y T ,Z T ) 

with Xjn be the total wealth at time T in bonds and index 5*. In turn, the value function 
reads 

V{t, x, y, z) = sup E \u {X T + U a {Y T , Z T )) X t = x,Y t = y, Z t = z] (3) 

where M. is a set of admissible trading strategies that require holding of initial cash x. 
The expectation in the Eq.([3]) is taken under the "real- world" measure P. 

For a portfolio made exclusively of stocks/indices and bonds, the value function for 
the exponential utility is known from the classical Merton's work: 

V°(x,t) = - e -^ xerT -^ T (4) 

where r = T—t, r is the risk free interest rate assumed to be constant, and rj s = {^ s —r)/a s 
is the stock Sharpe ratio. 

In our setting, in addition to bonds and stocks/indices, we want to long Cz option 
and short a units of Cy option to statically hedge our Cz position, or, equivalently, buy 
the g% option. 



The stock is equivalent to our index S in the setting of the Merton's optimal investment problem. 
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The value function in our problem of optimal investment in bonds, index and the 
composite option has the following form: 

V(x, y, z, t) = sup E f- e -^T+n«(Y T ,z T ) ^ = X}Yt = y, Z t = z) (5) 

where Xt is a cash equivalent of the total wealth in bonds and the index at time T. We 
represent it in a form similar to Eq.Q: 



1„2 



V(x,y,z,t) = -e-^ e -^$(y,z,T) (6) 

where function $ will be calculated in the next sections. The indifference pricing equation 
reads 

V{x,y,z,t) = V°{x + g%-ap Y ,t) 
Plugging this in Eq.Q and Eq.([6]) and re-arranging terms, we obtain 

9z = --e rT log$(y,z,T) + ap Y 

7 

The highest price of the Z-derivative is given by choosing the optimal static hedge given 
by the number a of the F-derivatives, i.e. 

9°z = ~-e rT log $ a . (y, z, r) + a*p Y (7) 
7 

a* = arg max < — e rr log z, r) + apy 
a I 7 

where we temporarily introduced subscript a in $ a to emphasize that the value function 
depends on a through a terminal condition. 



3 The HJB equation for HIMD 

To use Eq.([7]) and thus be able to compute both the option price and optimal static 
hedge, we need to find the "reduced" value function $. To this end, we first derive the 
Hamilton- Jacobi-Bellman (HJB) equation for our model, and then obtain its analytical 
(asymptotic) solution. 

Let 7r = 7Tf(x) be the dynamic investment strategy in the index St at time t starting 
with the initial cash x, and C n be the Markov generator of price dynamics corresponding 
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to strategy tt. Both the optimal dynamic strategy and the value function should be 
obtained solution of the HJB equation 

V t + sup £"V = (8) 

7T 

We assume that all state variables St, Y t , Z t follow a geometric Brownian motion pro- 
cess with constant drifts pi and volatilities (Ji,i G (x,y,z) 

dSt = p x Stdt + a x StdWt 
dY t = p y Y t dt + a y Y t dwi y) 
dZ t = p x Z t dt + a z Z t dW t (z) 



If our total wealth at time t is X t = x and we invest amount n of this wealth into the index 
and the rest in a risk-free bond, the stochastic differential equation for X t is obtained as 
follows: 

dX t = r(X t - tt) dt + ^dS t = (rX t + ira x r) s ) dt + na s dw} x) , rj s = ^— - 

St cr x 

Then C* reads 

C n = (rx + ir(n x - r)) V x + -a 2 7r 2 V xx + p y yV y + -crffiVyy + +p z zV z 
+ ^a 2 z z 2 V zz + p xy cr x cr y nyV xy + p xz cr x cr z TczV xz + p yz (T y a z yzV yz , 

where V(x, y, z, t) is defined on the domain WL(x, y, z, t) : [0, oo) x [0, oo) x [0, oo) x [0, T] . 
Since C n is a regular function of tt, sup n is achieved at 

* / \ f]s^x Pxii@i{yV xy ~\~ p xz o z zV xz 
it Jx) = " " 



@xV xx 



Plugging this into Eq.ph, we obtain 



V t + rxV x + p y yV y + -a y y V yy + p z zV z + -a z z V zz + p yz a y a z yzV yz (9) 



1 (VsVx + PxyV y yVxy + Pxz^ z zV xz ) 2 _ 

2 yir.fr. 
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This is a nonlinear PDE with respect to the dependent variable V(t, x, y, z) with standard 
boundary conditions (see [15]) and the terminal condition determined by a choice of the 
writer's maximal expected utility (value function) of the terminal wealth Wt- 

Note that so far the derivation is valid for a generic utility function. To make further 
progress we specialize to the case of exponential utility in Eq.(]2j since it gives rise to a 
natural dimension reduction of the HJB equation. Indeed, the ansatz 

V(t, x, y,z) = - exp (- 7 a;e r(T - T) ) G(h, s, r) (10) 

with h = log(y/K y ), s = log(z/K z ) is both consistent with terminal condition Eq.(j5| and, 
upon substitution in ([9|, leads to a PDE for function G which does not contain variable 



x: 



1 Jlr, , 1 _2 



G T — flyGh + £l z G s + -CyGhh + ^ z Gss + PyzVy&zGhs (11) 



2 

1 2(~i 1 {PxyVyGh + PxzCzGs) 2 

~ 2 Vs ~ 2 G ' 

Here 

1 2 . 1 2 

Py = Py- ^ ~ VsPxyOy , Pz = Pz~ Tfz ~ VspxzVz 



Equation Eq.(ll) is defined at the domain ~R(h,s,t) : [—00,00) x [—00,00) x [0,T]. The 
initial condition for this equation is obtained from Eq.Q. 

In what follows, we choose a specific payoff of the form Eq.([T]) with ITy = min(Y, K y ), Hz 
min(Z, K z ) that corresponds to a portfolio of bonds of firms Y and Z with notionals K y , K z 
within the Merton credit-equity model. Then the terminal condition for G(h, s, r) reads 

G{h, s, 0) = exp [-7 (K z e s - - aK y e h -)] (12) 

where s_ = min(s,0) and h_ = min(/i,0). 



4 Asymptotic solutions of Eq.(ll) 



We were not able to find a closed form solution of Eq.(ll) with the initial condition 
Eq.(12). On the other hand, a numerical solution of this equation is expensive, espe- 
cially when it should be used many times for calibration to market data. Therefore, we 



proceed with aasymptotic solutions of Eq.(ll). We suggest two approaches to construct 
asymptotic solutions. 
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4.1 First method 



As we want to statically hedge option Cz with options on another underlying, we look 
for an asset Y that is strongly correlated with asset Z. Further, if we have a "similar" 
option Cy on asset Y (i.e. similar maturitytype, strike, etc.), we expect that such option 
provides a good static hedge for our option Cz Q 

Therefore, a natural assumption would be to consider 1 — p yz to be a small parameter 
under our setup. Utilizing this idea we represent the solution of Eq.(ll) as a formal 
perturbative expansion in powers of e: 



(13) 



8=0 



where e is the small parameter to be precisely defined in the next section. 



As we shown below, Eq.(ll) can be solved analytically (in quadratures) to any order 



of this expansion, thus significantly reducing the computation time. 



4.1.1 The HJB equation in "adiabatic" variables 

We start with a change of variables (h, s) —>• (w, v) defined as follows: 



w 



a,, 



PxzO'z 



_i_ T ( —On. _|_ ^ z P xy 

@y "z Pxz 



Simultaneously, we change the dependent variable G — > $ as follows: 



G(h, s, r) = e 2^ T $(w, v, r) 



Using Eq.(|ll|), Eq.(|14|) and Eq.(]15|), we obtain the following PDE for function $: 

$2 



1 



i /4 

2' ww 2 



^PxyPxzPyz H" P X y 
rxz 



* 'j 



Pxz PxyPyz i 
T ^ t 

Px2 



1 



\Px 



XL) 



(14) 



(15) 



(16) 



Further we will show that v is a slow ("adiabatic") [^variable of our asymptotic method, 
while w becomes a "fast" variable. 



5 As an example, we mention the case of equity options referencing the same underlying, i.e. Y = Z, 
but K z K y . We may want to hedge an illiquid option with strike K z (say, deep OTM) with a liquid 
option on the same underlying but with a different strike K y . Under this setup, we have p yz — 1, i.e. a 
prefect correlation case. 

6 See e.g. [? ]. 
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In what follows we need an inverse of Eq.(14) at r = 0: 

Pxz<?z 



h = fjyiv , s = (v + w)- 



Pxy 



Using this in Eq.(12), we obtain the initial condition in (w,v) variables for the function 

v, r): 

$(tu, v, 0) = exp -7 (K z e%£ a * {v+w) - - aK y e a * w A (17) 
where (x)_ = min(x, 0) for any real x. 



4.1.2 Cosine law in 3D and "adiabatic" limit 

Recall that a correlation matrix £ of iV assets can be represented as a Gram matrix 
with matrix elements Ejj = (xj, x^) where Xj, Xj are unit vectors on a iV — 1 dimensional 
hyper-sphere S N ~ 1 . Using the 3D geometry, it is easy to establish the following cosine 
law for correlations between three assets: 



Pxy = Pyzpxz + \J {l - p 2 yz ) (1 ~ pL) COS(<f) xy ), 



;is) 



with (j) xy being an angle between x and its projection on the plane spanned by y, z. 

As discussed e.g. by [2], three variables p xy , p yz , (p xz are independent, but p xy , p yz , p xz 
are not. Therefore, one of them, e.g. p xy , has to be found using Eq. ( 18 ) given p xz , p yz , (p xy . 
Further we define e as 

(19) 



1-P^<1, 



and also define the following constants 

0i = 



I -P. 



xz 



COS 



Pa 



J xy)i 



1+01 



(20) 



Using Eq.(18) and definitions in Eq.(19), Eq.(20), coefficients at Q uv and <& vv in Eq.(16) 



are evaluated as follows: 

Pxy Pyz 



-1 + 



P, 



e @3i Pxy — PxzPi 



Plz ~ 2p X ypxzPyz + P. 



P 2 

Hxy 



p = vi -£ 2 + e6i, e 3 = Vi- e 2 e x - e. 



Accordingly using this notation Eq.(16) takes the form 

1 



$r = -<S> ww + e6z§ 



E 2 6 2 $ vv 



2 Hxzl $ 



(21) 



In the limit £ — > this equation does not contain any derivatives wrt v, therefore v enters 
the equation only as a parameter (since G(w,v,r) is a function of v). We call this limit 
the adiabatic limit in a sense that will be explained below. 

It should be noted that our expansion in powers of e can diverge if p xy is very small. We 
exclude such situations on the "financial" grounds assuming that all pair- wise correlations 
in the triplet (S t ,Y t ,Z t ) are reasonably high (of the order of 0.4 or higher in practice), 
for our hedging set-up to make sense in the first place. Thus, parameter 9\ is treated as 
0(lfl 



4.2 Second approach 



It turns out that the last equation of the previous section could be further simplified. 
Introducing new independent variables 



1 



u 



2 w 



0-i 



v/e 



(22) 



we can transform Eq.(21) into the following equation 



(23) 



It is seen that in new variables the mixed derivative drops from from the equation, as 
so does e. However, further let us formally introduce a multiplier \i in the term which 



transforms Eq.(23) into 



1 ^ ^ 



(24) 



Let us also formally assume that [i is small under certain conditions. The idea of this trick 



is as follows. One way to construct an asymptotic solution of the Eq.(23) is to assume 



that all derivatives are of the order 0(1), and then estimate all the coefficients. If one 
manages to find a coefficient which is 0(/x), then it is possible to build an asymptotic 
expansion using that coefficient (or fi itself) as a small parameter. If, however, all the 



While parameter 62 is always O(l) and positive, parameter 6\ could be both positive and negative 
for typical values of correlations. For example, if (p X y, pxz: Pyz) = (0.4,0.4,0.8), then Q\ = 0.33, while 
for (p xy , Pxz, Pyz) = (0.3,0.2,0.8) it is 0\ = —0.22. The cosine law can also be used to find proper 
values of correlation parameters in the limit p yz — >• 1. To this end, we first use Eq.(18l to convert the 
estimated triplet (p X y, Pxz-, Pyz) m to a triplet of independent variables (p xz , 4> X y, Pyz), and then take the 



limit 



Pyz 



1 while keeping p xz and cf) xy constant. 
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coefficients in the considered PDE are of order 0(1) we need to check if perhaps some of 
the derivatives in the Eq.(|23| are small, e.g. O(p). If this is the case, in order to apply 
standard asymptotic methods we formally have to add a small parameter p as a multiplier 
to the derivative which is 0(p) f\ make an asymptotic expansion on p, solve the obtained 
equations in every order on p, and at the end in the final solution put p = 1. That is 
exactly the way we want to proceed with. 

This means that instead of the Eq.((T3|), we now have the following expansion: 



(25) 



To find conditions when & m could be small as compared with the other terms in the 
Eq.(23), we use an inverse map at r = : (h, s) — > (u, v) 



(J,, 



03 



+ u 



+ u/3 



and rewrite the payoff function Eq.(12) in the form 

<$>(u,v,0) = exp [-7 (K z e Ca * {uJ1+u) - - 
9i 



0* 



UJ 2 



aK y e^ y ^ 2+u) -)] 

1 



(26) 



^2 /3v^2 V$2 

Suppose that v > 0, u < —uj\ or v < 0, u < —U2- Differentiating the payoff twice 
by u and twice by v and computing the ratio of the first and second terms in the rhs of 
the Eq.(23), one can see that in the limit e — > this ratio becomes \i = O^vv/^uu — ®\- 
Typical values p xz = 0.3, p xy = 0.2, p yz = 0.8 give rise to \i = 0.05, therefore the second 
term is small as compared with the first one, and p is a good small parameter. This, 
however changes if p xz is small or/and cos((p xy ) is close to 1, and then p oc 0(1). Still in 
this case we have e < 1 which can be used as a small parameter. Therefore, our approach 
is as follows: 



1. If Q\ 1 we use Eq.(23) and find its asymptotic solutions using Eq.(25). This is 
better than using Eq.(21) because first, p is typically smaller then e, and second, 
the term Q m i n our first method Eq.(13) appears only in the second order of ap- 



proximation while in the second method it is taken into account already in the first 
order on p. 



8 In other words write $ ss = fi(^vv / /i) 
9 At e — ?> this is a regular map. 



fi(^vv), where $ BS oc O(l) 
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2. If, however, Q\ oc 0(1), then p is not anymore a small parameter, therefore we use 



Eq.(21) and solve it asymptotically using Eq.(13) 



In general, this argument cannot be applied if v > 0, u > — w 2 or v < 0, u > —uj\ 
because then both derivatives of the payoff vanish. However, the above argument is 
intended to provide an intuition as to why Q m could be much smaller that the other 



terms in the rhs of Eq.(23). This intuition can be verified numerically, and our test 
examples clearly demonstrate that smallness of p often takes place. Below we discuss 
under which conditions this could occur. 

Note that at the first glance, the described method looks similar to the quasi-classical 
approximation in quantum mechanics ([13]). The similarity comes from the observation 



that transformation Eq.(22) is singular in e which is similar to the quasi-classical limit 



h — > 0. If we would construct an asymptotic expansion on e we would expand the rhs 



of the Eq.(23) on e, but not the payoff function. After getting the solution of Eq.(23) 
in zero-order approximation on e as a function of {u,v), we would apply the inverse 
transformation (u, v) — > (h, s) which is non-singular. Therefore, the final result would not 
contain any singularity. Since p = 0\ is defined via p xz , p xy and p yz = a/1 — e 2 , it could 
seem that p = p(e), and we face a "quasi-classical" situation. 

However, as explained above, the independent parameters are p yz ,p xz and cos((f) xy ). 
Therefore, by definition p = Ol(p xz ,cos((f) xy )) doesn't depend on p yz , or on e. Thus, our 
two methods actually correspond to different assumptions. The first one utilizes a strong 
correlation between assets Z and Y. The second assumes a strong correlation between 
index S and asset Z while at the same time the vector of correlation p xz in 3D space is 



not collinear to the vector of correlation p xy . By financial sense this means (see Eq.(18)) 
the following. 

1. Either p xz is about 1 and, therefore, p xy ~ p yz . In other words, index S strongly 
correlates to asset Z, so S is almost Z, therefore correlation of Y and Z (p yz ) is 
close to correlation of Y and X {p xy )- That, in turn, means that asset Z can be 
dynamically hedged with S, and extra static hedge with Y doesn't bring much 
value. In contrast, under the former assumption static hedge plays an essential role. 

2. Or cos((p xy ) <^ 1 which means that p xy ~ p xz p yz . For instance, p xz = 0.4, p yz = 
0.6, pxy — 0.24. This is an interesting case, since it differs from two previously 
considered assumptions on high value of either p yz or p xz . Indeed, all correlations 
could be relatively moderate while providing a smallness of Q\. 

In what follows, we describe in detail the asymptotic solutions for zero and first order 
approximations in p, and outline a generalization of our approach to an arbitrary order 
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in [l. Asymptotic solutions in e are constructed in a very similar way and are given in 
Appendix |Aj Also to make our notation lighter in the next section we will us v instead of 
v since that should not bring any confusion. 



4.3 Zero-order approximation 



In the zero order approximation we set /i = 0, so that Eq.(23) does not contain derivatives 
wrt v: 

* - l *> l tf ( $ 0^) 2 (27) 

where p\ z = 2 p 2 z . Therefore, dependence of the solution on v is determined by the 
terminal condition. In other words, our system changes along variable u, but it remains 
static (i.e. of the order of /x 2 slow) in variable v. Using analogy with physics, we call this 
limit the adiabatic limit. 

The last equation can be solved by a change of dependent variable (closely related 
to the Hopf-Cole transform, see e.g. Henderson & Hobson [5], Musiela & Zariphopoulou 

my- 

$ (r, u, v) = [<p(r, u, v)] 1/(1 - pL) (28) 



which reduces Eq.(27) to the heat equation 

1 



2 ruui 



(29) 

subject to the initial condition <f> U V Q = v, 

0y~p'L. The latter can be obtained from 
the Eq.(26) if one replaces 7 with the "correlation-adjusted" risk aversion parameter 
7 = 7(1 — p 2 2 ). It can also be written as a piece- wise analytical function having a 
different form in different intervals of w-variable. If v > 0, we have 

{exp [-7 {K z e^^ 1+u) - aK y e^ u y^ 2+u) )] , u < -u x 
exp [-7 (K z - aKye^y^ 2 ^)] , -wi < u < -w 2 (30) 

exp [-7 (K z - aKy)} , u > -u 2 , 

while for v < we have 

{exp [-7 [K z e Ca ^ 1+u) - aK y e^ a y^ 2+u) )] , u < -u 2 
exp [-7 (K^^ 1 ^ - aK v )] , -u 2 < u < -u t 

exp [-7 (K z - aK y )} , u > -u x 
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Using the well-known expression for the Green's function of our heat equation Gq(u' — 



u,r) 



(l£ —u) 



v see e.g. [IE]), the solution of Eq.(36) is then 

(f)(u,v,r) 



V2 



TIT J — 



e 2t <p{u , v, 0)du 



The explicit zero-order solution thus reads 

1 



$ (u,v,t) 



IV2 



TTT 



00 (,'-„)» ' 

e 2r (f)(u ,v,0)du 

-00 



V(i-pL) 



(31) 



Note that Eq.(31) provides the general zero-order solution for the HJB equation with 



arbitrary initial conditions at r = 0. For our specific initial conditions Eq.(30), the 



solution is readily obtained in closed form in terms of the error (or normal cdf) function 
(see Appendix [B]). However, for numerical efficiency it might be better to use another 
method which is based on a simple observations that the expression in square brackets 



in the Eq.(31) is just a Gauss transform of the payoff. This transform can be efficiently 



computed using a Fast Gauss Transform algorithm which in our case is 0(2N) with N 
being the number of grid points in u space. 



4.4 First-order approximation 



For the first correction in Eq.(23) we obtain the following PDE 

2 



$0,u 



uu Pxz ^ 



(32) 



with Qi(u,v,t) = ^2 < l ) o,™/2. This is an inhomogeneous linear PDE with variable coeffi- 



cients. As long as our zero-order solution of Eq.(27) already satisfies the initial condition 



this equation has to be solved subject to the zero initial condition. This considerably 
simplifies the further construction. 



We look for a solution to Eq. (32) in the form 



$i(m, v, r) = [$ (m, v, t)\ p " H(u, v, t) 

This gives rise to an inhomogeneous heat equation for function H subject to zero initial 
condition 



;H UU + 0!$ 



14 



Thus, using the Duhamel's principle Q16J) we obtain 

(u-U 1 ) 2 

<$>i{u,v,t) = $%"(u,v,t) / dx du' = 

Jo J-oo v / 27r(r-x) 



® x {v!,v,x)®7 x *{v!iV,x) (33) 



There exists a closed form approximation of the internal integral (see Appendix D). 



4.5 Second order approximation and higher orders 



The second order equation has the same form as the Eq.(32) 



O.u 



$ 2 + e 2 (M,w,r) 



where 



6 2 (u,t;,t) = ^6» 2 $i iTO - ^ 



-2 



$1$0,„ o $l$0,u$l,u , 

$o 



$3 



*2 



As this equation has to be solved also subject to zero initial conditions, the solution is 
obtained in the same way as above: 



$ 2 (u,v,t) = $%"(u',v,t) 



e 2 ( t -x) 



=e 2 (it, v, x)$ Pxz («'» v i X)dxdu' 



'o J-oo v/27r(r - x) 

This shows that in higher order approximations in \i both the type of the equation and 
boundary conditions stay the same. Therefore, the solution to the n-th order approxima- 
tion reads 



t roo 



e 2 ( t -x) 



=0„(tt, v, x)$ Pxz K> w ' X)dxdu' 



lo J-oo \/2ix(t - x) 

where G n can be expressed via already found solutions of order i,i = l...n — 1 and their 
derivatives on u and v. The exact representation for n follows combinatorial rules and 
reads ^ 

n (tl,V,r) = -02^n-l,vv ~ ^Plz-n, 

where S n is a coefficient at p n_1 , n > 1 in the expansion of 



L l + /3 2 (/i) 



O.u 

oo 



8=1 



15 



in series on //. This could be easily determined using any symbolic software, e.g. Mathe- 
matica. For instance S3 reads 

„ (-$l$0,n + $Q$l, u ) (gfgo^ - $Q$l$l,u + 2$p (-ggggg + $0$2,n)) 

The explicit representation of the solutions of an arbitrary order in quadratures is impor- 
tant because, per our definition of /i, convergence of Eq.(Jl3|) is expected to be relatively 
slow. Indeed, if one wants the final precision to be about O(0. 1) at p yz = 0.8 (/i ~ 0.36), 
the number of important terms m in expansion Eq.po"]) could be rawly calculated as 
fi m+1 = 0.1, which gives m = 1.25, while at precision 0.01 this yields m = 3.5. 

Note that all integrals with n > 1 do not admit a closed form representation and have 
to be computed numerically. Again, this could be done in an efficient manner using the 
Fast Gauss Transform. 



5 Validation of the method and some examples 

To verify quality of our asymptotic method we compare two sets of results. One is obtained 
using zero and first order approximations (being computed via a series representation and 



f2 functions given in Appendix B in Eq.(41), Eq.(42) and Eq.(45), or using the Fast 



Gaussian Transform). Our tests showed that the number of terms in the double sum that 
should be kept is small, namely truncating the upper limit in i from infinity to i ma x — 10 
produces nearly identical results. Therefore, the total complexity of calculation is about 
45 computations of exp and Erfc functions which is very fast. A typical time required for 
this at a standard PC with the CPU frequency 2.3 Ghz ranges from 0.68 sec (Test 1) to 
0.35 sec (Test 2) (see below). 



The other test is performed using a numerical solution of Eq.(23). In doing so we use 
an implicit finite difference scheme built in a spirit of p2]. After the original non-linear 
equation is discretized to obtain the value function at the next time level, we need to solve 
a 2D algebraic system of equations each of which contains a non-linear term. This could 
be done e.g. by applying a fixed point iterative method (|17j). In other words, at the first 
iteration as an initial guess we plug-in into the non-linear term the solution obtained at 
the previous level of time. This reduces the equation to a linear one since the non-linear 
term is explicitly approximated at this iteration. Next we solve the resulting 2D system of 
equations with a block-band matrix using a 2D LU factorization. At the second iteration, 
the solution obtained in such a way is substituted into the non-linear term again, so again 
it is approximated explicitly. Then the new system of linear equations is solved and the 
new approximation of the solution of the original non-linear equation is obtained. We 
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Test 


P>x 


0~ x 


r 


Pyz 


K g 


Pz 


<y z 


Pxz 


^0 


Ky 


0v 




Pxy 


?/o 


1 


0.04 


0.25 


0.02 


0.8 


110 


0.05 


0.2 


0.4 


100 


90 


0.03 


0.3 


0.3 


100 


2 


0.04 


0.25 


0.02 


0.8 


110 


0.05 


0.3 


0.3 


50 


90 


0.03 


0.3 


0.2 


100 



Table 1: Initial parameters used in test calculations. 



continue this process until it converges. The number of iterations to needed for numerical 
convergence depends on gradients of the value function, which are considerably influenced 
by the value of 7. For small values of 7 (about 0.03) we need about 1-2 iterations, while 
for 7 m 0.3, 5-7 iterations might be necessary. For higher values of 7, the fixed point 
iteration scheme could even diverge, so another method has to be used instead. It is also 
important to note that we solve the non-linear equation using the dependent variable 
log($o), rather than $ to reduce relative gradients of the solution. 

Note that for linear 2D parabolic equations with mixed derivatives more efficient 
splitting schemes exist, see e.g. [9]. In principle, such schemes could be adopted to solve 
Eq.(23). Though Eq.(23) is a non-linear equation, the presence of the non-linear term 
does not change its type (it is still a parabolic equation), and second, does not affect 
stability of the scheme if we approximate it implicitly and solve the resulting non-linear 
algebraic equations. A more detailed description of this modification of the splitting 
scheme of Hout and Welfert will be presented elsewhere. 

Implementation of the numerical algorithm is done similar to [lOj. In typical tests we 
use a non-uniform finite difference grid in u and v of size 50x50 nodes, where —50 < u < 
50, —30 < v < 30. The number of steps in time depends on maturity T, because we use 
a fixed time step St = 0.1 yrs. Typical computational time for T = 3 yrs on the same PC 
at 7 = 0.03 is 17 sec. In Fig[l]a 3D plot of the value function V(u, v) is presented for the 
initial parameters marked in Table [T] as 'Test 1'. We also use 7 = 0.03, a — 1. It is seen 
that V(u, v) quickly goes to constant outside of a narrow region around u = and v = 0. 



In Fig[2]the same quantities are computed for vq = 1.22 and T =10 yrs. Here the first 
plot presents comparison of the numerical solution with zero and '0+1' approximations. 
The second plot compares the zero and first order approximation. It is seen that the first 
approximation makes a small correction to the zero one in the region closed to w = 0. 
Also both '0" and '0+1' approximations fit the exact numerical solution relatively well. 
This proves that our asymptotic closed form solution is robust. Results obtained with a 
second set of parameters (Test 2 in the Table [T]) are shown in Fig [3] In Fig [1] we present 
price g% computed in tests 1,2 as a function of u, where py is Black-Scholes put price 
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with parameters of the corresponding tests. 

Note that for values of parameters used above the nonlinear term is small, therefore 
the solution is closed to the solution of the linear 2D heat equation obtained from the 
Eq.(21) by omitting the nonlinear term. That is because 7 is small in our Test 1. To 
make testing more interesting, we changed 7 to 7 = 0.2. Results obtained for T = 3 yrs 
are presented in Fig [5] for Test 1. and in Fig [6] for Test 2. It is seen that in this case the 
'0+1' analytical approximation still fits the numerical solution. 

The computational time in these tests is higher because the matrix root solver con- 
verges slower. The typical time at T=3yrs and a = 1 is 24 sees. Computation of the 
analytical approximation requires the same time as before which is about 1 sec. 

If the exponent of the payoff function is positive, e.g. at a = 2 and other parameters 
as in Test 1, then the solution looks like a delta function. Under such conditions any 
numerical method experiences a problem being unable to resolve very high gradients 
within just few nodes. Therefore, in this case one has to exploit a non- uniform grid 
saturated close to the peak of the value function. This brings extra complexity to the 
numerical scheme while our analytical approximation is free of such problems. 

On the other hand it is natural to statically hedge Cz option with the option Cy with 
same or close moneyness. This significantly reduces the exponent and allows one to use 
the same mesh even at high values of the risk aversion parameter 7. 

As the numerical performance of the model depends on the value of 7, a few comments 
are due at this point. Though we do not address calibration of our model in the present 
paper (leaving it for future work), the value of 7 should be found by calibrating the model 
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Figure 2: Value function V{u,v) at T=10 yrs obtained by FD scheme, '0' and '0+1' 
approximations. The initial parameters are given in Table [T| test 1. 




Figure 3: Value function V(u,v) at T—W yrs, comparison of the zero and first order 
approximations. The initial parameters are given in Table [T| test 1. 
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Figure 4: Price g% obtained in tests 1,2. 




Figure 5: Value function V(u, v) and price g% at T=3 yrs, comparison of the numerical 
solution, '0' and '0+1' approximations at 7 = 0.2. The initial parameters are given in 
Table [II test 1. 
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Figure 6: Value function V(u, v) and price at T=3 yrs, comparison of the numerical 
solution, '0' and '0+1' approximations at 7 = 0.2. The initial parameters are given in 
Table [L| test 1. 



to market data. It is not entirely obvious how to do this since we use an illiquid asset Z, 
and moreover the price of our complex option g% is not an additive sum of its components 
for incomplete markets. It therefore makes sense to calibrate the model to another set of 
instruments that are both liquid and strongly correlated with the original instruments. In 
the context of equity options, such calibration of 7 was done in p], giving rise to values 
of 7 ~ 0.1 — 0.6. While it is not exactly clear how 7 found in this setting is related to 7 
of our original problem, we expect the latter to be of the same order of magnitude. 

Note that the results in Fig|6] clearly demonstrate that the proposed method is just 
asymptotic. Indeed, the first correction to the solution in Figj6] is a few times larger 
then the zero-order solution. As usual, there exists an optimal number of terms in the 
asymptotic series that fit the exact solution best. We do not pursue such analysis in this 
papei^ 

To characterize sharpness of the peak in the value function it is further convenient 
to introduce a new parameter ip = ^K z . If ip is high, e.g. tp > 50 the value function is 
almost a delta function, so the asymptotic solution as well as its numerical counterpart 
are not expected to produce correct results unless they are further modified. Based on 
the results of [1] one can see that the value of ip calibrated to the market varies from 
0.01 to 15. Therefore, in our test 2 we used ip = 20. Both our numerical and asymptotic 

10 From theory we know that the asymptotic solution converges to the true solution as e — > 0, however 
due to our definition of e via p yz this is not an interesting case. 
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methods work with no problem for these values of ip. 

In Fig [7] the optimal hedge en* is computed based on Eq. (j7|) which was solved using 
Brent's method ([E]). The initial parameters correspond to test 1 in Table [l] in the first 
plot, and test 2 in the second plot. Note that for u < —15 (the first plot) and u < —20 (the 
second plot), Eq.(j7]) does not have a minimum, so the maximum is obtained at the edge 
of the chosen interval of a. The latter could be defined based on some other preferences 
of the trader, for instance, the total capital she wants to invest into this strategy etc. 

Finally, in Fig [8] price g§ is presented as a function of a for various 7. It is seen that 
this function is convex which was first showed in [8] in a different setting. Note that 
these results were obtained using a new numerical method mentioned above. It combines 
Strang's splitting with the Fast Gaussian Transform, and accelerates calculations approx- 
imately by factor 40 as compared with a non-linear version of the 2d Crank-Nicholson 
scheme. A detailed description of the method will be given elsewhere. 
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Figure 8: Price function of a at various 7 computed for the initial data in tests 1 

and 2. 
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A Asymptotic solutions of the Eq.(21) 



Here we describe in more detail the solutions for the zero and first order approximation 
in s, and outline a generalization of our approach to an arbitrary order in e. 



A.l Zero-order approximation 



In the zero order approximation e = the Eq.(21) does not contain derivatives wrt v: 

\2 



$0,r = -^0, U u - 77/4 



(*o,«)' 







(34) 



The solution of this equation proceeds along similar lines to Sect. 4.3 using a change of 
dependent variable 



o(T,u,v) 



IT. 



u,v)} 



V(i-rt.) 



which reduces Eq. ( 34 ) to the heat equation 



subject to the initial condition 4> U)V> q = v, 0) 



(35) 



(36) 



The explicit form for the latter 



coincides with (30) provided we substitute 7 = 7(1 — pl z ), co± = 1, £ = l/{3 and U2 = 0. 



The explicit zero-order solution thus reads (compare with Eq.(31)) 



$ (u,v,t) 



V2 



ITT 



~4>(u', v, 0)du 



V(i-pL) 



(37) 



Note that Eq.(37) provides the general zero-order solution for the HJB equation with 



arbitrary initial conditions at r = 0. For our specific initial conditions Eq.(30), the 



solution is readily obtained in closed form in terms of the error (or normal cdf) function 
(see Appendix [Bl . 



A. 2 First-order approximation 



For the first correction in the Eq.(21) we obtain the following PDE 

2 



where Qi(u, v,t) = 9 3 & 



<I.. -p 2 <1> "- !i >'< ' 1 



XZ Q 







(3f 
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This equation coincides with Eq.(32 ) except that the free term is different, and the cor- 



relation parameter is p xz rather than p xz . The solution proceeds as in Sect. |4.4[ resulting 
in the following expression: 



(v 



Q 1 (u,v,t) = Q*>(u,v,t) I d x I dv! ; 2(T X) =Q x {v! ,v,x)®/ xz {u! ,v,x) (39) 

J0 J-oo y/2n{T-X) 

The double integral that enters this expression can be split out in two parts. One of them 
could be found in closed form, while the other one requires numerical computation (see 
Appendix [C]) . Our numerical tests show that for many sets of the initial parameters the 
first integral is much higher than the second one, so the latter could be neglected. However, 
we were not able to identify in advance at which particular values of the parameters 
this could be done. Moreover, for some other initial parameters we observe an opposite 
situation. 



A. 3 Second order approximation and higher orders 



The second order equation has the same form as Eq.(38) 



where 



e 2 (u, v, t) = -9 2 ®a, m + #3 - #3$o,™) - ^pL$o ($i/$o)L 2 



As this equation has to be solved also subject to zero initial conditions, the solution is 
obtained in the same way as above: 



$ 2 (u,v,t) = V xz (u',v,t) / / --Q 2 (u,v,x)® Pxz (u ,v,x)dxdu' 

J0 J-oo V 27r ( r ~ X) 

This shows that in higher order approximations in e the type of the equation to solve 
doesn't change as well as the initial conditions. Therefore, the solution to the n-th order 
approximation reads 



2 r e 2( T - x ) _ 2 

$ n (w, v, r) = « v,r) / Q n (u, v, x)$ Pxz v, xWdu 1 

Jo J-oo V 27r ( T - X) 
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where G n can be expressed via already found solutions of order i,i = l...n — 1 and their 
derivatives on u and v. The exact representation for n follows combinatorial rules and 
reads 



Q n (u,v,r) = 3 ^2 



1 







tn — i+l 



i=i 



n 



Dl Q e n-i+l 



i—l,uv 



(40) 



e=0 







n-i+2(2 



i=2 



n-i + 2)\ de n ~ i+2 



i—2,vv 



e=0 



where S n is a coefficient at e n , n > 1 in the following expansion 



<91n$ c 

du 



;i + /3 3 



i + 



dln(l + /3) /dln$ 



(9m 



du 



-i 



/3 3 = $>^M 



i=l 



In particular, S 3 reads 



—3 



+ 



0,U 



+ 2$ 2 , u 



The explicit representation of the solutions of an arbitrary order in quadratures is im- 
portant because per our definition of e the convergence of the Eq.(13) is expected to be 
slow. Indeed, if one wants the final precision to be about 0.1 at p yz = 0.8, the number 
of important terms m in the expansion Eq.(13) could be rawly calculated as e m = 0.1, 
which gives m = 4.5. 

All the integrals with n > 1 do not admit a closed form representation and have to be 
computed numerically. 



B Closed form solutions for the zero-order approxi- 
mation 



Since our payoff is a piece-wise function, the integral in Eq.(31) can be represented as a 
sum of three integrals. We denote u 
follows: 



0{U,V,T) 



4« + jp 



pxyv/pxz and represent the zero-order solution as 
C = sign(w) 



J. 



(0 



27 



where sign (— ) means that u = p xy v/p xz < 0, and sign (+) - that uj > 0, and 



(+) 



J. 



{-) 



V2 



ITT 



du e ^ exp 



-7 [K z e »*y - aK y e y J 



^ i-u,-n/K z e zp *y ,cr z — ,>yaK y ,a y \, 

V Pxy J 

1 f° , (u'-^f r / a 
y du e 2r exp -7 (^ - aK y e° yU J 

= Q (0, -7X2, 0, 70;^, (7^) - Q (-co, -jK z , 0, jaK y , a y ) , 



,--y(K z -aK y ) roo 



V2 



1 - -Erfc 
2 



72 



7TT 7- 



dw'e 2t exp 



( 0, -7^/ J *»V z — ,~/aK y ,a y ) , 



1 



y/2 



ttt Jo 



_ (u -u) 

du e it exp 



^[K z e ff '^ u+u,) -aK. 



= D, ( -w, -^K^'^ ,a t — ,7Q;if y ,0 ) - D, ( 0, -^K z e" ^ w , a z — , ~/aK y , ) , 



(41) 



V2 



TTT 



Pxy 

du'e-^ = V^-^Erfc 
2 



■u + 00 



Here Erfc(x) is the complementary error function, and 



Q(a,5,p,p,q) 



V2 



g e pu',gqu' _ (u'-u) 2 
3 +P e 2r 



7TT J-oo 

00 i 



1 ?TlfP_ 



M i - j)+qj]u 'e- (J! ^du' 



i-Z7XjKi-3) } - V v/27 / 



i=o j=o J v J7 
A = ip + j(g — p). 



(42) 



Since the complementary error function quickly approaches zero with 1 < 0, or 1 with 
x 3> 0, the number of terms one has to keep in the above sums should not be high. 
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We also need the following function 



1 f°° (u'-uf 

Ju,v = r- — / e 2r <p u , )V (u', v, 0)du' , (43) 



T 



OG 



when computing the first order approximation. Using a general form of the payoff 
(f)(u, v, 0) = e s ( v ) ePU +P equ [iij it can be represented in the form 



J u ^ = -j== / e 2- (f> u > !V (u',v,0)du' (44) 



IDC 

OO 



where 



<r»^L= / e~ i ^e p "'0(w>,O) fp + e 9U 'g/3 + e^'p^u)) du' 

V27TT J-oo V / 

T (C) , T (C) , T (C) 



_ 7-C+) _ t-(-) _ n 

<^2,u,v <^3,u,v <^3,u,v u ' 

T(+) - 7^ I -rf Vz^U Pxz - ts \ 

Ji,u'v = ~lK z a z e **v iii I —lo, —jK z e <>*■» ,a z — ,jaK y ,a y \ 

V Pxy J 

t(-) - ts ^f^n fn -rs vz^u Pxz - is \ 

Ji,u'v = -lK z a z e r*y \li I 0, —^K z e »*y ,a z — , jaK y , a y I , 

V Pxy J 

Jtul = -lK z a z e a ^\ fix f-u;,-7i^ e ^V*— ,7^,0 

V Pxy 



( 0,-7^ 2 e ffz ^V 2 ^,7a^,0 

Pxy 



and 



OiCa, tf.p, 0, q) = ^ ■u- J f- )l iP A (P) + P 6A ( 2 P) + ?0 A (P + ?)] > ( 45 ) 
i=o j=o 3>' 

A(6) = e^K^fcfc ^-^W^ ) s ip + j{q _ p) + 6 



Compare with Eq.(30) 
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C Transformation of the first order solution of the 



Eq.(21) 



The first order approximation is given by Eq.(33) where the zero-order solution <&o(u, v, r 



has been already computed in Appendix B We plug in this solution into Eq.(33 ) to obtain 



$i{u,v,t) =$g"(u, v,t) / d X du' =e 1 (u',v,x)®o Pxz (u',v,x) (46) 

Jo J-oo a/27t(t — x) 



J 



2 1 r~ e i^-x) - p *z 1 

6 3 ® p *°(u,v,t) / d X / riu' /n . - J{u\v,x) TZ &d <v J(u',v,x) TZ % 
Jo J-oo y/2n(T-x) 



where the integrals J, 



(0 



1, 3 are defined in Eq.(41) 



The internal integral can be simplified. Indeed, since 



i-p: 



2 



:i-pD 2 j 



the internal integral in Eq.(46) can be represented as a sum of two integrals where 

(u — it ) 

1 f°° . , e"^ 3 ^ 



i-p: 



£2 •/ — 00 



P 2 

rxz 



2 ^2 

xz) 



1 r Ju'vi 

VMr-x) 

(u-u') 2 

, e 2 (-x) J^J M , 
2vr(r - x ) ^ 



As we already mentioned the second integral could be either smaller or larger than 
the second one depending on the parameters. This is illustrated in Fig. [9] The initial 
parameters are taken from test 1 in Table [l] with a = 1 and T = 3 yrs. In Fig. [10] the 
same calculation is shown for the parameters corresponding to test 2 in Table [TJ 
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Figure 9: Comparison of 3\ and 3 2 at T=3 yr and 7 = 0.03 and 0.2, test 1. 



1st approximation - T-3 years, v - -2.87S994e+00Q 
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Figure 10: Comparison of 3i and 3 2 at T=3 yr and 7 = 0.03 and 0.2, test 2. 
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The first integral 3% can be rewritten using the integral representation of J in Eq.(31 ) 



1 f°° e~ 2 ^-^ f e ^ — 

du' d u i I du" — 4>o :V (u", v, 0) 



1 - Plz i-oo ^2tt(t - X ) U J-oo V^X 



(u-u') 2 (u"-u') 2 

1 f°° , „ ± . „ . f°° , , e^^^T u "-u'e ~ 



4>o,v(u ,v,0) I du 



1 - Plz J-oo ' ' ' 7-oc ^2vr(r - x) X 
1 f°° e 

- 1 - f i ii i /ii ^ \ ^ 



Substituting this into Eq.(46) and integrating over x, we find 



$ 1 (u,v,t) = 1 -^-$$"(u,v,t) d X du"^ v {u\v^)-=—{u" -u) (47) 



1-/4 U Vo ' ' 

~, C«-«") 2 

= — 4-^«(tt,v,r) / du"<j) 0tV (u", v, 0)^-=^(u" - it) 

= - r 9sT T ^ L (u,v,r)J u>v . (48) 

Thus, the first correction to the solution obtained in the zero order approximation on e 

is approximately — 9^t \_ 2 yz &n xz (u, v , t)J UV) and the full solution in the "0+1" approx- 
imation reads 



v, r) = $ (u, v, r) - ggT ^ (u, v, r)J UjV 

1 ~Pi 



xz 



D Transformation of the first order solution of Eq. ( 23 ) 



The first order approximation is given by Eq.(33). It is assumed that the zero-order solu- 



tion $o( M ; v , T ) is already computed by using either the method presented in Appendix |Bj 
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or using the Fast Gauss Transform. We plug in this solution into Eq.(33) to obtain 



$i(u,t;,t) = $ p * z (u,v,t) I d X 
1 



du 6 2(T ' X) . 81K, v, X )% Plz («', v, X) (49) 



g 2(t-x) -Prz 1 

dw' — =J(u', v, x) 1 ~^ z d V)V J(u' } v, x) l ~ p * z 



J 



J® + Jf 



J. 



(0 



where the integrals J, 



(0 



1,3 are defined in Eq.(41) 



It can be seen that Eq.(49) is similar to Eq.(- 



6) if one replaces p xz with p xz and 



d VjV J(u', v, x) 1 ~ p '* z with d VjV J(v, v, x) 1 ~ p ^ g ■ Therefore, we can use the same idea as in 
Appendix O to further simplify this integral. 



Accordingly the inner integral in Eq.(49) can be rewritten as 

1 



1-A 



Jy.v ~4 



plz 



o ... , 



2 ^2 



J 



The internal integral in Eq.(49) can then be represented as a sum of two integrals 3i + 3 2 , 
where 



1 



1 - P 2 , 

rxz J — oo 
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P 2 

rxz 



Get . J V ^ VJ 

v 27r ( r - X) 

(u, — u ) 



27r(r-x) J 



The first integral Ji can be modified using the integral representation of J in Eq.(31) 
1 



3i 



1 - p 2 . 

x rxz J — oo 



00 , . e~^^> 
du 



du"- 



■2\ 



1 " Pi 
1 



VMr ~ X) 
du"(f)o,vv(u" ,v,0) I du - 



xz J — oo 



^x 

(u-u') 2 

e 2 ( r - x) e 



0o,w(w",i;,O) 



2X 



'27r(r - x) 



Pi 



du"(p o ,vv(u",v,0)- 



V2 



TIT 
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Substituting this into Eq.(49) and integrating over x, we find 



$i(u,v,t) = - 2 Z#'(u,v,t) / d X du"(j) 0>vv (u",v,o f (50) 

1 1 - Pxz J0 J-oo \J2lXT 

OO (u-u") 2 

-$q xz (u,v,t) / du"(f> o>vv (u",v,0)- 



2 1 - J-oo V27TT 

Thus, the first correction to the solution obtained in the zero order approximation on 

-2 

/i is approximately -tOi tz^z &q xz (u, v, t) J v v , and the full solution in the "0+1" approx- 

J- Pxz 1 

imation reads 

$(u,v,t) = $ (u,v,t) + ^rfl 2 _^_ 2 $f*(u,v,T)J v , v . 

^ L Pxz 
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